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In this paper we present a comparison of trajectory optimization approaches, 
for the minimum fuel rendezvous problem. Both indirect and direct meth- 
ods are compared for a variety of test cases.. The indirect approach is 
based on primer vector theory. The direct approaches are implemented 
numerically and include Sequential Quadratic Programming (SQP). Quasi- 
Newton, and Nelder-Mea.de Simplex. Several cost function parameteriza- 
tions are considered for the direct approach. We choose one direct approach 
that appears to be the most flexible. Both the direct and indirect methods 
are applied to a variety of test cases which are chosen to demonstrate the 
performance of each method in different flight regimes. The first test case 
is a simple circular-to-circular coplanar rendezvous. The second test case is 
an elliptic-to-elliptic line of apsides rotation. The final test case is an orbit 
phasing maneuver sequence in a highly elliptic orbit. For each test case we 
present a comparison of the performance of all methods we consider in this 
paper. 

INTRODUCTION 

The minimum fuel rendezvous problem has received considerable attention in the literature and 
numerous approaches to both posing and solving the problem have been developed. The basic objective 
is to find a minimum solution to a two-point boundary value problem (TPBVP) that has multiple feasible 
solutions. The primary difference between approaches is the choice of independent variables and how 
the TPBVP is posed in terms of these variables. Methods also differ on how constraints are handled. 
In the next few paragraphs we give a brief overview of the methods considered in this paper. We also 
discuss the test problems we used to allow a comparison of each method we consider for different flight 
regimes. 

This paper deals with the problem of optimizing trajectories where the maneuvers can be modeled 
impulsively. In general, two approaches are utilized for optimization: direct and indirect . 1 Both ap- 
proaches will be compared in this investigation. A survey of the rendezvous problem is provided by 
Jezewski et al . 2 

A direct optimization method tries to minimize the cost function by directly varying the control ( in - , 
dependent) variables. The methods employed for the optimization process are mathematical program- 
ming techniques . 3 The mathematical programming techniques utilized in this paper include Sequential 
Quadratic Programming (SQP), Quasi-Newton, and Nelder-Mea.de Simplex. (Other techniques, such as 
Genetic Algorithms, are under current investigation but are not discussed here). For direct optimization 
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the choice of independent variables and constraints is extremely important- For example, suppose that 
there are n impulses in the trajectory of interest. Then, two possible choices for independent variables 
are (1) the actual AV's, or (2) the impulse locations and times. In this paper, we find the option of using 
the impulse locations and times more effective. This is consistent with previous work by Brusch. 4 Hence, 
we utilize a Lambert-type approach with implicit position continuity. Thus, as the impulse locations 
change, new arcs are computed utilizing either Lambert arcs (in the two-body model) or differentially 
corrected arcs (allows for more complex models than the two-body model). The result of computing the 
lambert arcs is a path continuous trajectory with velocity discontinuities. The velocity discontinuities 
are then utilized to formulate the cost function. Other approaches are possible but include the addition 
of position continuity as constraints. 5 

An indirect optimization method tries to minimize the cost function by considering its variations, 
thus it involves Calculus of Variations ( COV). It leads to two-point boundary value problems in the 
co-state variables . Once the co-state variables a, re obtained, the optimal control variables are usually 
readily available. Primer vector theory can be considered to be a byproduct of applying the Calculus 
of Variations to the problem of minimizing the fuel usage of impulsive trajectories. Perhaps more 
importantly, primer vector theory can be utilized to “find' 5 the optimal number of impulses. A survey 
of primer vector theory is presented by Guzman 6 et ai 

The parameterization of the cost function chosen in this work ensures that for every cost function 
evaluation, the rendezvous constraints are satisfied implicitly. Hence, from the point-of- the- view of a 
numerical optimizer, the problem is unconstrained. Several numerical optimization routines are con- 
sidered. Three of the routines are products of The Mathworks and are available in their Optimization 
Toolbox. The Mathworks’ routines that we employ are fmincon, fminunc, and fminsearch. The fmincon 
function is an SQP algorithm, fminunc is a Quasi-Newton algorithm, and fminsearch is a Nelder- Meade 
Simplex algorithm. For a detailed discussion of The Mathworks 5 routines and their implementations, 
we refer the reader to The Mathworks Optimization Toolbox 7 documentation. We also use SNOPT, an 
SQP algorithm developed by Gill 8 et. ai For details on SNOPT we refer the reader to the user’s guide 
for SNOPT 8 5.3. For the indirect method we develop software in Matlab that uses a fully automated 
algorithm which implements the primer vector theory first developed by Lawden. 9 

We employ three test cases to compare the performance of each trajectory optimization method 
in different flight regimes. The first test case is a simple circular- to-circular, coplanar rendezvous. The 
second test case is an elliptic-to-elliptic, line of apsides rotation. The final test case is a phasing maneuver 
in a highly elliptic orbit. For each test case we generate numerous initial guesses. The initial guesses are 
generated using simple intuition. The times and locations of the initial and final maneuvers are chosen 
based on experience. Given the initial and final burn locations we generate a two-burn rendezvous 
sequence. Multiple maneuver sequences are generated by placing small maneuvers, equally spaced in 
time, along two maneuver initial guesses. Sequences of two, three, and four burns are considered. The 
performance of all methods for the test cases described above is presented in the final section. It is 
difficult to provide a single metric that can conclude which method is the best. Hence we have provided 
a variety of statistics to illustrate the performance of all of the trajectory optimization approaches we 
consider. It should be noted that we have not considered the rate of convergence to be a measure 
of importance. Since all methods are acceptably fast, we use metrics based on the total AV of the 
converged solution to be the primary measures of performance. 

To solve the minimum fuel rendezvous problem we must first develop a mathematical model. In the 
next section we discuss the distinction between an orbit transfer and an orbit rendezvous and develop a. 
model for the rendezvous problem. 

PROBLEM STATEMENT 

The minimum fuel rendezvous problem can be posed in numerous ways. The primary difference 
between approaches is the choice of independent variables, how the boundary value problem is posed in 
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Terms of the independent variables, and how the constraint functions are handled. It is useful at this 
point, to make a clear distinction between a rendezvous problem and a transfer problem. In an orbit 
transfer problem we are only concerned with finding a maneuver sequence that will take a spacecraft 
from some position in its initial orbit, to some position in the desired final orbit. For a transfer, the 
initial maneuver can occur a ny where along the initial orb i t . an d t h e hi i a 1 m a 1 1 e u ve r can occ ur an y w h ere 
along the final orbit, without concern for the orbit phasing. However, for a rendezvous problem, we have 
the additional constraints determined by the time and place of the spacecraft in the initial orbit, and 
the time and place of the target spacecraft in its orbit. In this section we develop a mathematical model 
of the rendezvous problem. In later sections the model is used to develop optimization algorithms to 
find minimum fuel solutions. 

A diagram to aid in illustrating the general impulsive rendezvous problem is shown in Figure i. The 
arc denoted V 0 represents the path of the spacecraft in its initial orbit. Similarly, the arc denoted V f 
represents the path of path of the spacecraft in its final orbit. Throughout the paper a superscript ' — ' 
sign is used to denote a condition immediately prior to an impulsive maneuver, and a superscript 
sign is used to denote a condition immediately following an impulsive maneuver. According to these 
conventions, we can write the conditions for the i th burn as, 

if = rf (1) 

f = (2) 

AV< = Vf-V+ (3) 

where rf and r~ are the positions immediately before and after the i th burn respectively, tf and i~ 
are the times immediately before and after the i th burn respectively, and V t + and V” are the velocity 
vectors immediately before and after the i th burn respectively. The constraints in Eqs.(l) and (2) are a 
result of the definition of an impulsive maneuver. 



Figure 1: Impulsive Orbit Transfer Diagram 


There are additional constraints that must be satisfied at the boundaries, and at the interior impulses 
which are imposed by the orbit dynamics. We assume that the dynamics for the entire maneuver sequence 
are given by 


r 



(4) 


Let Tf = f(r 0 . y 0 ) t 0 ,tf) and vj = g(r 0l v 0 , t 0 ,tf) be the solution to Eq. (4), with the initial conditions 
(r o ,v 0 ,£ o ) and the final time of £/. For an interior impulse i ^ 1, i n, and n > 2. The additional 
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constraints that must be satisfied at an interior impulse are given by 


ry = f(r y.j.y) (5) 

v t + = ( 6 ) 

The constraints that must be satisfied at the initial boundary are given by 

r i" = f(r 0 ,v 0 ,t. 0 J. l ) (7) 

v t = g(r o-.Vo-1o-ti) (8) 


where (r 0 . v 0 ,£ 0 ) are the position and velocity of the initial orbit at the reference epoch t Q , and t\ is the 
time of the first maneuver. A similar set of constraints apply at the final boundary and are give by 

(9) 

( 10 ) 

where (ry,vy,fy) axe the position and velocity of the initial orbit at the reference epoch tf, and t n is 
the time of the last maneuver. Without loss of generality* we can appropriately formulate the optimal 
control problem to ensure that the constraints r* = r~ and tf = tj are satisfied implicitly. For the 
remainder of this work we drop the superscripts on r and f. and assume 

r * = = r r 

U = tt= r- 

In summary, the problem is to find the set (r^, xf . v~, £*), where i = 1.2...n and n > 2, that satisfies 
the constraints 


ri 

= f(r<>,v 0 ,t 0 ,*i) 

(13) 

v i" 

= g(r o ,v 0 ,t 0 ,ti) 

(14) 

ri 

= f(r t _i : v i _ 1 ,f i _ 1 ,t i ) 

(15) 

v t r 

= g(r»-i 

(16) 

r n 

= f(r f ,v f ,tf,t n ) 

(17) 

v+ 

v n 

= g 

(18) 


and that is a minimum solution to the function 

j=y;iiAv i n (19) 

i=i 

where n is the number of maneuvers, the set (r 0 ,v 0 ,£ 0 ) defines the initial orbit and the set (r /, vy,£y) 
defines the final orbit. 

There are two main methods available for solving the optimal control problem, including direct and 
indirects approaches. Within each method there are numerous possible implementations. Often methods 
are actually a combination of indirect and direct techniques. For direct methods, the parameterization of 
the problem depends on the type of optimizers available. Indirect methods are sensitive to the constraints 
and often require complete reformulation if new constraints are imposed. In the next two sections we 
discuss several approaches for solving the minimum fuel rendezvous problem using both indirect and 
direct methods. 

DIRECT APPROACH 

The success of a direct method is intimately dependent on the choice of independent variables and 
the constraint function formulation. In this section we discuss three parameterizations for solving the 


(11) 

( 12 ) 


v- = g(r f.vj-.tf.t,,) 
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minimum fuel rendezvous problem using direct methods. There are numerous possible parameterization? 
and comparing all approaches is beyond the scope of this work. We categorize the parameterizations 
we consider into two groups called the Feasible Iterate Approach and the Infeasible Iterate Approach. 
We present a brief comparison of the paiameterizations. and the most promising method is chosen for 
further development. For this work we choose a Feasible Iterate Approach for reasons explained later. 
This section is concluded with a discussion of several numerical optimization packages used in this work. 

Objective Function Parameterizations 

There are numerous approaches to solving the minimum fuel rendezvous problem using direct meth- 
ods. For convenience we break down the direct approaches considered here into two categories called the 
Feasible Iterate Approach and the infeasible Iterate Approach. In the Feasible Iterate Approach, each 
evaluation of the objective function is also a feasible solution that satisfies the rendezvous conditions. 
Hence, from the point of view of the optimizer, the problem is unconstrained. The Feasible Iterate 
Approach requires that for each objective function evaluation TBPVP is solved. The Feasible Iterate 
Approach also requires a careful selection of independent variables for the objective function parameteri- 
zation. In the Infeasible Iterate Approach, there is less restriction on the choice of independent variables, 
however, the constraints must, be defined carefully, and the optimizer must be able to handle nonlinear 
constraints. For continuity, we present the discussion of alternative parameterizations in Appendix 1. 
Here we present the parameterization chosen for this work. 

Our Choice of independent variables is as follows: 

Given: (r 0 ,v o ,t 0 ) and (r y,vy,i/) (20) 

Choose the independent variables: 

t-i i = 1,2.. .n ( 21 ) 

Tj j = 2, 3....n — 1 (22) 

where t* are the times of the maneuvers, and rj are the positions of the intermediate burns. Given the 
independent variables defined in Eqs.(21) and Eqs.(22), we must develop an algorithm to determine the 
total AV r of the maneuver sequence. The boundary conditions ri, vj% r n , and v~ are determined using 
Eqs.(13),(14),(17),and (18), where ti and t n are known from Eq.(21). After solving for the boundary 
conditions we know the times and positions of all the maneuvers. Therefore, we can solve for the 
velocities before and after each maneuver, by solving Lambert’s problem for each of the n — 1 trajectory 
segments. There are numerous well-known approaches to solving Lambert’s problem. x We choose a 
method developed by Howell and Peraicka 10 and further developed by Guzman. 11 The method developed 
by Howell is chosen because it can solve Lambert's problem with perturbations included, as well as when 
there is more than one significant gravitational body. Perhaps more importantly, Howell demonstrated 
that using the method described below, we can provide analytic approximations for the gradient of the 
total AV with respect to the independent variables chosen in Eqs. (21) and Eqs.(22). The algorithm is 
described as follows. Given an initial position r*_i, and an initial velocity both at time U-i, find 
applied at time so that we achieve r * at U* Figure 2 illustrates the problem. The dark black 
arc denotes the path a spacecraft would follow if no £v*_i is applied. For this arc, the final position 
denoted by r a is not equal to the desired final position r t . Hence £r t ^ 0. The dashed arc denotes the 
trajectory that is the solution to Lambert's problem. For this arc r; = r 0s or St{ = 0. To solve for 
such that 8ti = 0 first define x as 



then <5x t is given by 

6xi = $ (tuti-i)6xi-i (24) 
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where t, arid £,•_ i are fixed. We can write # (r t , — r ) as 


Using the fact that 5rj_i — 0 we can write 




<5ri = 


Solving for i we see that 
From inspection of Fig. 2 we can write 
Solving for Sr f we have 
Substituting Eq.(29) into Eq.(27) yields 


<5v,_! = B“ _ l( 5ri 
r* + Jr, = r„ 
Sri = r a - r< 


Jvi-! = B { ,_1 (r a - r<) 


(26) 

(27) 

(28) 

(29) 

(30) 


We can solve for such that Sr l is zero by iterating on Eq.(30) until (r 0 — r, ) meets a user defined 
tolerance. It is important to note that this approach assumes that B^_ 1 exists. For cases when B~ i _ 1 
is not invertible we use a method by Gedeon 12 to solve the TPBVP. We solve all n - 1 trajectory arcs 
using the algorithm defined above, then the AV is calculated using Eq.(19). Note that although the 
algorithm above can be extended to provide a gradient approximation, we have used finite differencing 
for gradients in this work. Analytic approximations for the gradient are. a topic of current research. 


There are several other concerns to address to completely define the objective function parameteriza- 
tion. We also must choose appropriate time and coordinate systems to express the independent variables 
given in Eqs. (21) and (22). There are several issues in choosing a time parameterization that must be 
considered. The first issue is to select appropriate units, the second is to select an appropriate reference 
epoch. We choose the units of seconds for time, and reference each time to the reference epoch to 
given in Eq. (20). The positions in Eq. (22) are expressed in cartesian coordinates in the Mean of J2000 
Earth Equatorial system and the units are in km. 


With the parameterization of the cost function described above, the rendezvous conditions are sat- 
isfied implicitly for every cost function evaluation. Hence, the problem is essentially an unconstrained 
optimization problem and we are free to use unconstrained as well as constrained optimization techniques 
to find a minimum AV solution. It is exceedingly rare that both TPBVP methods described above fail to 
converge. In the case that neither method converges we must “inform” the optimizer. This is discussed 
in a subsequent section. In the next subsection we discuss the numerical optimization packages we use 
in this work. 
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Direct Optimizers 


There are numerous optimization packages available that solve unconstrained optimization problems. 
We investigate the performance of four routines for the parameterization of the minimum fuel rendezvous 
problem chosen above. Three of the routines are products of The Math works and are available in their 
Optimization Toolbox. 7 The first routine is an SQP algorithm called fmincon. The second routine is a 
quasi-Newton algorithm called fminunc. The third routine is a Nelder-Meade simplex algorithm and is 
called fminsearch. For a detailed discussion of the optimization routines developed by The Math works 
we refer the reader to the Optimization Toolbox documentation/ The fourth numerical optimization 
routine we employ is an SQP routine developed by Gill 8 ei al. called SNOPT. For a detailed discussion 
of SNOPT we refer the reader to the SNOFT 8 5.3 documents t km. 


Numerical Issues 


To avoid numerical difficulties the independent, variables and the objective function values are normal- 
ized to be on the order of one. Secondly, all derivatives are calculated using finite differencing. Providing 
analytic derivatives for gradient based methods is a topic of current research. The final numerical issue 
occurs when both TPBVP solvers fail to converge. This is exceedingly rare. However when it occurs we 
must “inform” the optimizer. SNOPT has a built in capability to step back from a poorly conditioned 
state vector. Hence, for SNOPT, if both TPBVP solvers fail, the appropriate message is sent to SNOPT 
and the SNOPT steps away from the poorly conditioned state.’ However, The Mathworks routines do 
not have this capability yet. As a temporary solution, when using a Mathworks’ routine, in the event of 
a poorly conditioned state vector we pass back a crude approximation of the total AV. 

To solve the minimum fuel rendezvous problem using direct methods we must choose an appropriate 
parameterization of the problem, and an appropriate numerical optimization routine. In this section 
we discussed some choices involved in choosing a specific parameterization. We chose one method that 
performed well in a preliminary comparison. We also discussed several numerical optimization routines 
we employ and some numerical issues one must deal with to ensure adequate performance. In the next 
section we discuss an indirect approach to solve the minimum fuel rendezvous problem. 


INDIRECT APPROACH 

In this section, we present a review of the primer vector theory as well as the main challenges 
involved in its implementation. 13 Primer vector theory 14 is an optimization technique based on calculus 
of variations. The theory has several appealing features including the indication of when and where to 
add an impulse to a non-optimal trajectory, and a visual assessment of the optimality of neighboring 
solutions. 


Primer Vector Theory 


Primer vector theory provides a set of first order necessary conditions, which a trajector}' must meet 
to be locally optimal. The necessary conditions, first derived by Lawden, are expressed in terms of 
the primer vector, which is defined as the adjoint to the velocity vector in the variational Hamiltonian 
formulation. 14 If any of Lawden’s conditions are violated, the rendezvous trajectory is not optimal and 
we can use the primer vector history to obtain information on how to improve its AV r cost. In his initial 
work, Lawden solved a fixed-time rendezvous and his theory was further extended by Lion and Handels- 
man and la, ter, by Jezewski and Rozendaal to solve the n-impulse optimal rendezvous problem. 15 ’ 16 A 
more detailed derivation of the primer vector theory is provided by Hida.y. 17 For the rendezvous problem 
considered in this paper, we can express the primer vector equations using calculus of variations theory. 
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The initial trajectory or first-guess, is labeled as the reference Trajectory. Since primer vector is a first- 
order theory based on local variations, it will converge to local optimal neighboring trajectories of the 
reference trajectory. Therefore, the optimal solution is highly dependent on the reference trajectory bur 
it will also depend on other design parameters discussed later in this section. The primer vector obeys 
the following equation, also known as the second order canonical form of the Euler- Lagrange equation. 

p=/n-^+ 3 ^(I>vU) (3i) 

i=i 

where r is the satellite position vector on the reference trajectory, p is the primer vector and p is the 
Earth's gravitational constant. The satellite position is found using Eq. (4). As evident in Eq. (31). 
the primer vector state car not be derived simultaneously with the spacecraft state. The spacecraft 
trajectory must be propagated first for the primer vector history to be computed. Using calculus of 
variations. Law den derived four necessary conditions expressed in terms of the primer vector defined by 
Eq. (31) for an optimal rendezvous trajectory: 


1. The primer vector must be C 1 (i.e. the primer vector and its first derivative are continuous) for 
the entire history. 

2. Between impulses. |jp|j = p <— 1. Impulses occur when p = 1. 

3. At an impulse, p = p = u^. where u? r is the optimal thrust direction. 

4. At all interior impulses (not at the initial or final times) p • p = 0. This condition has implications 
on the slope of the primer vector magnitude since 


<*HpII 

dt 




P P 

iipir 


Therefore, = 0 at the intermediate impulses. Also, for convenience, let p — 


These conditions are necessary (but not sufficient) for an optimal trajectory. In this paper, a trajectory 
that meets the above conditions will be called an optimal trajectory. When solving a rendezvous problem 
using primer vector theory, we first need to evaluate the primer vector history along the reference 
trajectory or initial guess. Solving for the primer vector history is equivalent to solving a TPBVP. Let’s 
first assume that we are solving a two- burn rendezvous problem. We know from Lawdems conditions 
that, at an impulse, the primer vector is in the direction of the thrust vector. Thus, we can express the 
initial and final primer vector as : 

p 0 = T^rr. (32) 


p/ = 


I A Vo] ’ 
A V/ 

|AVf| 


(33) 


The initial and final primer vectors and the time-of-flight (£/ — t 0 ) are known. However, to propagate 
Eq. (31), we need the complete primer initial state (Po.Po). To obtain p 0 , we can either use a shooting 
method (TPBVP solver) or use the STM formulation from Eq. (34) below. 



A B 
C D 



(34) 


where A, B, C and D are 3x3 matrices, partitions of 3>(t,t 0 ). In general, using the STM is faster and 
although not as accurate as the TPBVP solver, it is sufficient. Losing Eq. (34), the initial derivative of 
the .primer vector can be expressed as: 


Po = B 1 • ( P / - A • p 0 ) 


(35) 
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It. is worthwhile noting that if B -1 is singular. p 0 cannot be computed using Eq. (35). In the two-body 
problem, there are three known singularities. 18 Eq. (35) can be generalized for a n-impulse trajectory 
by examining each individual two- burn transfer. Once the primer vector history is computed, we can 
determine the optimality of the trajectory using Law den's necessary conditions. If any of Lawden's 
conditions are violated, the rendezvous trajectory is not optimal and we can use the primer vecioi his lory 
to obtain information on how to improve its AV cost. Improving non-optimal trajectories using primer 
vector theory is the main contribution of Lion and Handelsman. 15 Essentially, first order variations of 
the total AV" are considered for different perturbed trajectories. 15 ' 17 For a non-optimal primer vector 
history, two types of actions are possible to lower the AV" cost: (I) moving the initial or final impulse and 
(II) adding and/or moving an interior impulse. If the epoch of the endpoints is unconstrained, an action 
of type I is performed whenever the initial and/or the final primer vector magnitude slope is different 
born zero. The general expression for time variations at both endpoints is given by Hi day 1 1 as 

SJ = - p 0 \\ Aw 0 \\dt o - j\f HAvyjjcff f, (36) 

where the initial and final time variations are given by dt 0 and dtj respectively. Furthermore, the initial 
and final primer magnitude slopes are given by p Q and pf, respectively. To have a lower cost, SJ must 
be less than zero. If the following terminology is defined, 1) initial coast (dt 0 > 0), 2) early departure 
(■ di 0 < 0), 3) final coast (dtf < 0 ). and, 4) late arrival (dtf >0 ): then, the following four cases cover 
all the non-zero primer slope combinations. 


• If p 0 > 0 and pf < 0 => Apply Initial Coast and Final Coast. 

• If p 0 > 0 and pf > 0 => Apply Initial Coast and Late Arrival. 

• If p 0 < 0 and pf < 0 =£ Apply Early Departure and Final Coast. 

• if p 0 < 0 and pf > 0 => Apply Early Departure and Late Arrival. 

See Hiday 17 for details. Of course, if the slopes are zero, no further improvement is achieved by varying 
the endpoints times. If the primer vector magnitude goes above one during a coasting phase on a given 
segment, an impulse is added to the trajectory (action type II) to improve the overall AV budget- This 
impulse is examined by considering a neighboring path with an additional impulse. The two-impulse arc 
is perturbed with the addition of a, mid-impulse at some time, t m , and position r m + <5r m . (The position 
r m is the position along the unperturbed path at t m .) The initial and the final times and positions, 
(t Q , r 0 ) and (tf, 17 ), are fixed. The cost function variation is defined as SJ = J’ - J, where J' is the 
fuel cost on the perturbed trajectory. Then, considering first order terms only, it can be shown that 

= C (1 — Pm*7)> ( 37 ) 

where the following are mid-impulse parameters: c is the magnitude of the impulse. p m is the primer 
vector at t m and 77 is an unit vector in the direction of the impulse. Then, for an improvement in 
cost, 5J < 0. Thus, using the definition of a dot product, if ||p m || > 1 at any time, a third impulse 
is beneficial. Furthermore, the greatest decrease in the cost function will be achieved if the impulse is 
applied at the maximum of ||p m |i at time t m and. in the direction of f). The position along the perturbed 
path and the magnitude of the impulse are yet to be determined. Consider finding the position of the 
impulse. Utilizing the fact that the position across an impulse is continuous, using the STM before and 
after the impulse, and utilizing the information obtained from Eq. (37), gives the following variational 
equation, 

*- c *' , iror 1381 

where 

A = D(t m , tf) 

— D(f m T 0 )B(^ m , 1 . 
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This equation, of course, is valid only for non-singular A and determines the position of the mid-impulse. 
To estimate the magnitude, c. of the impulse, the expression for the cost on the perturbed path. J', can 
be expressed as a function of c and minimized. The mid- impulse should decrease the cost but might not 
produce an optimal trajectory in the sense of Lawden. The subsequent optimization of the three- impulse 
trajectory is presented next. For a three-impulse arc. if thp time and position of the mid-impulse are 
allowed to vary, it can be shown that. 

6J = (Ap m) T dr m + (A H m ) dt m (39) 

where Ap rn and A H m are. respectively, the primer vector and the Hamiltonian differences at the im- 
pulse. If the trajectory is indeed optimal, the first cost variation must vanish. Thus, it is required that. 
Ap m = 0 and. A H m — 0. It can be shown that A H m — p~ T v ~ - p£ s T v+. In an interesting “mix" . a 
direct optimization method can be applied to vary the mid-impulse time ancl position to meet, the above 
conditions. In this investigation.- following Hidav, 1 ' the Broyden-F letch er-Goldfarb-Sh anno. (BFGS), 
algorithm is used. 


Primer Vector Implementation 


The* primer vector theory described above was implemented in Matlab in tool we call PVAT. PVAT has 
a fully automated algorithm, which iterates following the primer vector principles to optimize a non- 
optimaf reference trajectory. Most of the time multiple actions for improvements are possible. However, 
there is yet no mathematical theory that determines the optimal sequencing of the actions and, in gen- 
eral, different sequencing will lead to different neighboring solutions. One way of solving this problem is 
to combined Eqs. (36) and (39) to form a unique cost function gradient. This gradient is then used in the 
BFGS optimization algorithm to simultaneously move the endpoints times and the midcourse impulses 
position and time. A separate check for the addition of an impulse is then required. In this paper, we 
choose to implement a sequential algorithm where the endpoint times and the midcourse impulses are 
varied with separate PVAT iterations. Each time the underlying trajectory is changed, the algorithm 
recomputes the primer vector data and re-iterates until no improvement can be made. A conceptual 
flowchart of PVAT is shown in Fig. 3. 

First, the initial and final endpoint gradients are evaluated to check whether their values are greater 
than zero within some specified tolerance tt . Note that only one endpoint is moved at a time for a given 
iteration and the priority is given to the boundary with the highest gradient value. To determine the de- 
parture/arrival sta.te which corresponds to a zero initial slope, we implemente a bissection method. Both 
endpoint algorithms being identical, only the initial boundary is discussed here. Whenever the primer 
magnitude slope, is positive the departure state is propagated forward, otherwise it is propagated 
backward. To quickly recompute the trajectory data, we use a Lambert solver. The primer magnitude 
slope is then updated for the new trajectory and the process is repeated until the value of the slope is be- 
low the tolerance. Whenever there are internal impulses, the norm of the internal gradient is evaluated. 
The midcourse burn(s) are moved whenever the norm exceeds a specified tolerance e m . To move the 
midcourse impulses states, we use the BFGS minimization technique. This variant of the Quasi-Newton 
method uses an approximation of the Hessian matrix instead of its direct evaluation. The approximation 
formula is called a BFGS update. Finally, if the primer vector is above some maximum value, p threshold > 
during a coasting arc, an impulse is added to the trajectory. Note that theoretically the threshold value 
is one but for numerical implementation, the value is typically chosen to be slightly higher than one. The 
value of perturbation impulse is computed using an estimate derived by Jezewski. 16 However, we need 
to ensure that its magnitude remains small for the theory to be valid. In this paper, we implemented 
two different algorithms labeled PVAT1 and PVAT2 respectively. For PVAT1, the algorithm was not 
permitted to add a burn in spite of indication of potential improvements by the primer vector history. 
The second implementation permitts the addition of burns to the rendezvous trajectory to take full 
advantage of the primer vector theory. However, for practical reasons, we limited the number of burns 
to a maximum of six. For this first version of the code, the equations where not non-dimensionalized 
which means that the choices of the various tolerance are specific to the problem solved and that the 
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Figure 3: PVAT Conceptual Flowchart 


final optimal solution will depend on their value. 


TEST CASES 

In order to compare the performance of each optimization method we utilize three test cases. The 
test cases are chosen to investigate the performance of the each method in different flight regimes. The 
first test case is a low-Earth, circular- tcncircuiar, coplanar rendezvous. The second test case is an elliptic- 
to-elliptic line of apsides rotation. The third case is a highly elliptic orbit phasing sequence. For each 
test case there are numerous initial guesses. To generate different initial guesses for a particular test case 
we first define an initial and final orbit. Next, we vary the true anomaly and the epoch for the initial 
and final maneuvers. Given the times and positions of the initial and final maneuvers for a particular 
initial guess, we solve Lambert’s problem to yield a two burn solution. A simple algorithm is used to 
place small maneuvers along the two-burn trajectory arc to generate multiple maneuver solutions. In 
the next three subsections we discuss the three test cases in more detail. 


Case One 


Test Case One is a simple circular-to-circular coplanar transfer. The optimal solution is simply the 
Hohmann AV. The orbital elements for this test case are shown in Eqs.(40) and (41). The format for 
the orbital elements is [ a e % \ u H u ) unless specified otherwise. Distance is measured in km . , and angles 
are measure in degrees. 

Initial Orbit: oe 0 = [7000 0000 0], T = 0 (40) 

Final Orbit: oe f = [7500 0 0 0 0 M f ) , T —T f (41) 
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The AV for a Hohmann Transfer for Case One is 255.8 rn/s. A set of initial guesses for case one is 
generated by varying M f and Tf. We vary Mf from 90° to 270 c ' in increments of 20°. For each value 
of Mf we choose two values for Tj. For each pair of Mf and Tf we solve Lambert's problem to. obtain 
a two-maneuver rendezvous for the initial to the final orbit. The values of Tj are chosen so that the 
AV for the initial guess is less than 1.5 km/ s. Note this is considerably higher than the known optimal 
solution. For each two burn solution we add small interior burns to create similar three and four burn 
solutions. In summary, there are nine different values for A//, two different values of Tf for each Mf . 
and a two, three and four burn solution for each pair of Mf and Tf. Hence there are fifty four initial 
guesses for Case One. A figure showing a sample of the initial guesses for Case One is shown in Figure 

(4). 




Figure 4: Sample of Initial Guesses for Case One and Two 


Case Two 


The second test case involves a combined maneuver that changes both the semimajor axis and the 
inclination of the orbit. The orbital elements for Case Two are shown in Eqs. (42) and (43). The 
eccentricity of both orbits is 0.3. The semimajor axis for the initial orbit is 7000 km, and the semimajor 
axis for the final orbit is 8000 km. The initial orbit is equatorial, and has an argument of periapsis of 0°. 
Hence the five degree plane change is a rotation about the line of apsides. We generate 180 sets of initial 
conditions by varying the position and epoch of the initial maneuver, v 0 and T 0 respectively, as well as 
the position and epoch of the final maneuver, vf and Tf. We vary v Q from 300° to 40° in increments 
of 20°. We vary i/f from 140° to 220° in increments of 20°. These ranges are chosen because we know 
the initial maneuver should occur near periapis, and the final maneuver should occur near apoapsis. For 
each possible pair of v Q and Vf we choose two times of flight and compute two transfer trajectories that 
result in a two-burn maneuver sequence with a total AV of less than 2 km/s. Three and four burn 
maneuver sequences are generated from each two burn solution by putting small maneuvers, equally 
spaced in time along the initial two burn sequence. Hence, there axe thirty combinations of v 0 and Vf , 
two times of flight for each pair of v Q and i//, and two, three and four burn sequence for each time of 
flight resulting in one hundred and eighty initial guesses. A representative sample of the initial guesses 
for Case Two are shown in Figure (4). Positions of the initial maneuvers are label by black asterisks. 
Positions of the final maneuvers are marked by red circles. 

Initial Orbit: oe 0 = [7000 0.3 0 0 0 v 0 ] , T = T 0 (42) 

Final Orbit: oe f = [8000 0.3 5.0 0 0 u f ] , T = T> (43) 
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Case Three 


Case three is a phasing maneuver sequence in a highly elliptic orbit. The orbital elements for the 
initial and final orbits are given in Eqs.(44) and (45) where E 0 and Ef are the eccentric anomalies of 
the initial orbit and final orbit at times T 0 and Tj respectively. To generate a set of initial guesses for 
Case Three we vary E a from 0° to 345° in increments of 3 5°. For each value of E 0 we calculate two 
values of Ef such that there is a 30 minute separation in time between the initial and final orbits. The 
first value of Ef, for a given E a . is calculated so that the final position is thirty minutes ahead of the 
initial position. The second value of Ef . for a given E 0 , is calculated so that the final position is thirty 
minutes behind the initial position. For each Ej and E 0 pair we generate a two maneuver sequence by 
solving Lambert's problem for an Ef and E 0 pair knowing that At is thirty minutes. Three and four 
burn cases are generated by inserting small maneuvers along the two 'burn sequence. For Case Three 
there is a total of 96 initial guesses. 

Initial Orbit: oe 0 = [42000 0.8 10 0.0 0.0 E 0 ] , T = T 0 (44) 

Final Orbit: oej = [42000 0.8 10 0.0 0.0 Ef], T = Tf (45) 


The test cases developed above are intended to allow a comparison of the methods for rendezvous 
sequences in different flight regimes. In the next section we present results that compare the converged 
solutions for the six methods investigated in this paper. 


RESULTS 

Comparing the performance of the trajectory optimization methods studied in this work is non trivial. 
To begin, we must develop some metrics that allow. us to compare the solutions provided by the different 
techniques. The ultimate goal is to provide an analyst with some insight into which technique to employ 
for a given problem. So ultimately, we want to know which technique is most likely to converge to the 
lowest AV cost. However, this is not the only important statistic. We are also interested in determining 
the relative performance of different methods. For example, if one approach converges to a known 
optimal 70% of the time, and a second approach converges to a known optimal 25% of the time, the 
second approach may still be useful for several reasons. Firstly, for a given problem the hypothetical 
method 1 will not converge to the known optimal 30% of the time. Hence we will need another approach. 
Secondly, the hypothetical method 2 might converge to a slightly higher solution in terms of the AV 
cost, but it might provide a solution that is better when other mission constraints are taken into account. 

We employ two statistics to capture the relative performance of the different trajectory optimization 
techniques. The first statistic is simply the percentage of initial guesses where a particular method yields 
the lowest AV cost. We denote this statistic S\. For example, there are 54 initial guesses for Case One. 
If fmincon converges to the lowest AV cost of all the techniques, for 10 of the initial guesses, then Si for 
fmincon is 17.54%. According to this definition, the sum of the Si statistic for for all of the techniques 
for a particular test case is 100%. 

The second statistic we employ is intended to provide more insight into the relative performance of 
the different techniques. For example, if SNOPT converges to a solution of 260 m/s for a particular 
initial guess, and PVAT1 converges to a solution of 261 m/s for the same initial guess, then although 
PVAT1 has a higher solution, it still performed well relative to SNOPT. This is not captured in the S\ 
statistic. We define a second statistic S 2 , and define it as the percentage of initial guesses that converge 
to within 5% of the lowest solution yielded by the six techniques for the particular initial guess. An 
example makes this clear. Suppose we are comparing solutions from fminunc and fminsearch for a set of 
three initial guesses. Suppose fminunc converges to 100 m/s for the first initial guess, and fminsearch 
converges to 102 m/s for the first initial guess. For the second initial guess fminunc converges to 104 
m/s and fminsearch converges to 120 m/s . For the third initial guess fminunc and fminsearch converge 
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to 110. and i 0 1 m j s reaper t i vely. For this se i of solu t ion s 5 o it > r u n i nsea r cl i is 0 6 .00 % . which means tl \ a 1 
G G . G G % of the time f m insear eh conve rged t o within the 5 % of 1 1 1 e 1 owes t solution of ail the methods. 

Both the Sj and 5 2 statistic can be applied to the results in several ways. For example, if only 
four- burn initial guesses are used in generating an S -2 statistic, we add an additional subscript and call 
the statistic So 4 . If all of the initial guesses are used for determining an statistic, there is no additional 
subscript. This convention is also used for the S\ statistic. 

In the next three subsections we present a performance comparison of the six techniques for each 
of the three test cases. We use both the statistics described above, as well as graphical methods to 
illust rate how the techniques compare with one another. 


Case One 


Recall that Case One is a simple circular- to-circular copianar rendezvous problem. Hence the optimal 
solution is simply the Hohmann AY, which for the orbits defined in Eqs. (40) and (41) is 255.8 m/s. 
A figure illustrating the results for all six optimization techniques we investigate is shown in Fig. 5. 
The top subplot contains the solutions for all two burn maneuver sequences. The middle and bottom 
subplots contain the solutions for all three and four burn solutions respectively. On each subplot the 
x-axis is the initial guess number. The initial guesses have been numbered in order of increasing AV. 
This is done to enable one to draw conclusions as to the performance of each method as the AV of 
initial guess increases. Furthermore, the initial guesses between the subplots are similar. For example, 
the initial guess number one in the two-maneuver subplot is similar to initial guess number one in the 
three-maneuver subplot. The difference being only that a small mid-course maneuver is placed along the 
two-maneuver sequence to create a three maneuver sequence. The same is true for initial guess one of the 
four-maneuver subplot. On ail subplots the y-axis is the AV in km/s . Solutions for all direct methods 
are plotted in black. Solutions for indirect methods are plotted in red. Statistics that aid in comparing 
the results for each of method are found in Tablet. By inspection of Fig. 5 we see that the direct methods 
outperformed the indirect methods. For all but one initial guess, the direct methods always converged 
to within a few m/s of the known optimal solution for two-maneuver sequences. For three and four burn 
sequences, fminsearch consistently had the poorest performance of the direct methods. In comparisons 
between PVAT1 and PVAT2, PVAT2 performed better. Examining the S 1 statistics for Case One we 
see SNOPT found the lowest solution most often at 35.09% of the time. The method with the second 
best performance was fmincon, which found the lowest solution 31.58 % of the time. Examining the £2 
statistics we see that both SQP methods converged to within 5% of the lowest solutioin at least 93% of 
the time. It is also important to note that although other methods performed poorly in comparison to 
SQP, about 33% of the time the lowest solution was found by a method other than SQP. This suggests 
that although SPQ is the single best performing method for Case One, it is not safe to rely only on SQP 
as an optimization approach. It is also important to note that although SPQ found the lowest solution 
most often, solutions from other methods were often close to the solution found by SQP methods. 


Table 1: Statistics: Case One Results 



<Si 

Si, 

Si, 

•Si. 

S 2 

•S 2s 

•523 

<$2 4 

fmincon 

31.58 

21.05 

15.79 

57.89 

92.98 

100 

78.95 

100 

SNOPT 

35.09 

57.89 

42.11 

5.263 

94.74 

100 

89.47 

94.74 

fminunc 

10.53 

5.263 

5.263 

21.05 

78.95 

94.74 

52.63 

89.47 

fminsearch 

7.018 

10.53 

0.00 

10.53 

43.86 

100 

21.05 

10.53 

PVATl 

3.509 

0.00 

5.263 

5.263 

71.93 

94.74 

42.11 

78.95 

PVAT2 

12.28 

5.263 ' 

31.58 

0.00 

66.67 

57.89 

57.89 

84.21 
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Figure 5: Results for Case One 
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Case Two 


The maneuver sequence in Case Two involves raising the semi-major axis and rotating about the line 
of apsides, as well as satisfying the rendezvous conditions. Hence, two burn solutions will often perform 
poorly. Figure 6 shows all of the converged solutions for Case Two. The results axe plotted using 
conventions identical to the results for Case One and we refer the reader to the previous subsection 
for details. Table 2 summarizes the selected statistics for each method for Case Two. For the two- 
burn sequences seen in the first subplot, PVAT2 yielded the lowest AV about 48% of the time and 
converged 71% of time within 5% of the lowest solution. PVAT2 higher performance can be explained 
by its unique ability compared to the other methods to add maneuvers to lower the AV cost. For the 
two- burn sequences all methods, with the exception of PVAT2, frequently converged to the same local 
minimum. Note that PVAT1 did not perform as well as the direct methods for the two-burn sequences 
as it converged only 29% of the time within 5% of the best solution. For three-burn sequences, SNOPT 
and fmincon yielded to the lowest AV about 42% of the time and converged about 90% and 81% of the 
time, respectively, within 5% of the lowest solution. For this set of initial guesses, PVAT2 yielded the 
worst results. Finally, for the four-burn sequence, fmincon out-performed all the other methods, fmincon 
found the lowest solution 58% of the time and 93% of its converged solutions were within 5% of the lowest 
AV cost, fminunc is the second best with about 16% and 64% for Si and S 2 respectively. SNOPT, 
fminsearch and PVAT1 have comparable statistics for this initial guess sequence. PVAT2 performed 
poorly compared to the other methods. Overall, we can see that, for Case Two, when there are more 
than two-impulses, fmincon appears to be consistently better. However, for the two-burn sequence, 
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PVAT2 dominates over the other method because of its ability to add impulses to the trajectory. 
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Table 2: Statistics Case Two Results 



•Si 

<Si 2 

Si, 

Si. 

<s 2 

<52, 

<523 

<52, 

fmincon 

31.58 

9.677 

41.94 

58.06 

75.81 

53.23 

80.65 

93.55 

SNOPT 

24.73 

22.58 

41.94 

9.677 

56.99 

48.39 

90.32 

32.26 

fminunc 

6.452 

1.613 

1.613 

16.13 

50.00 

48.39 

37.1 

64.52 

fminsearch 

8.065 

8.065 

8.065 

8.065 

39.25 

43.55 

45.16 

29.03 

PVATl 

7.527 

9.677 

4.839 

8.065 

34.41 

29.03 

41.94 

32.26 

PVAT2 

16.67 

48.39 

1.613 

0.00 

33.33 

70.97 

17.74 

11.29 


Case Three 


Recall that Case Three involves a phasing maneuver in a highly eccentric orbit which implies that 
the initial guess is going to be critical in defining which optimal solution is achievable. This sequence 
is a pure phasing maneuver where the number of revolutions is not to exceed one. Because of the 
high eccentricity of the orbit, moving the initial and final epochs of the rendezvous will in some cases 
change the trajectory dramatically making it more difficult to set the proper tolerances for the finite 


16 




Total A V, km/s Total A V, km/s 


differencing to compute accurate gradients. Figure G shows all of the converged solutions for Case Three. 
Table 3 summarizes the selected statistics for each method for Case Three. For both the three- burn and 
four-burn sequences seen in the first subplot, all the indirect methods performed very poorly, fmincon 
outperformed all other methods with about 69%. and 66.7% of the lowest solutions for the three-burn 
sequence and the four-burn sequence respectively. For the three and four burn sequences. So for fmincon 
was 89.6%' and 9i.67% respectively, fminsearch and SNOPT were the direct methods with the lowest 
performance rate for both burn sequences. Note that fminunc converged 62.5% of the time within 5% of 
the best solution but was only the best method 10.42% of the time for the three- burn case. However, for 
the four-burn sequence, fminunc statistics improved to a convergence rate of 91.67% within 5% of the 
best solution and it. found the lowest AV 22.92% of the time. Overall, fmincon and fininunc seem to be 
the preferable methods for solving a Case Three type of rendezvous. All the other methods performed 
very poorly. 


TotalA V vs. Initial' Guess Number 
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Figure 7: Results for Case One 


Table 3: Statistics: Case Three Results 



Si 


Su 

s 2 

^2 3 


fmincon 

67.71 

68.75 

66.67 

90.63 

89.58 

91.67 

SNOPT 

11.46 

16.67 

6.25 

25.00 

25.00 

25.00 

fminunc 

16.67 

10.42 

22.92 

75.00 

62.50 

87.50 

fminsearch 

4.167 

4.167 

4.167 

9.375 

6.25 

12.50 

PVAT1 

0.00 

0.00 

0.00 

10.42 

8.33 

12.50 

PVAT2 

0.00 

0.00 

0.00 

10.42 

8.33 

12.50 
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CONCLUSIONS AND FUTURE WORK 


The minimum fuel rendezvous problem has received extensive attention in the literature. There are 
two primary methods to solving the problem: direct and indirect. However, many methods are actually 
a hybrid between a direct and indirect approach. 

In this work we compared several approaches for finding optimal rendezvous solutions. We compared 
several parameterizations for an objective function for a direct approach, and chose the one that appeared 
the best after some preliminary comparisons. We applied four numerical optimization routines to the 
direct method parameterization. It is important to note that finite differencing was used to calculate 
gradients for the gradient- based direct methods. Providing analytic gradients is a topic of current 
research and will likely improve both the rate of convergence, and the final values of the final converged 
solutions. 

We also developed two indirect approaches. The indirect approaches were based on primer vector 
theory. The first approach does not allow additional impulses higher than the number of impulses 
contained in the initial guess. The second approach will include an additional impulse if the optimality 
conditions suggest that it will result in a AV reduction. 

Ideally, from a software maintenance point of view, it is desirable for one method to always excel 
over competing methods. However, in this work we found that no single method was always the best. 
Yet, there are certain methods that outperform rival methods a significant amount of time. In general, 
in comparing direct approaches, SQP outperformed other methods such as Quasi-Newton, and Nelder- 
Meade Simplex. However, these results are intimately dependent, on the parameterization of the objective 
function. For alternative parameterizations, SQP might not necessarily outperform other methods. 

In comparisons between the direct and indirect approaches we implemented, we found that the 
direct approaches tended to find lower solutions more often. However, there were a significant number 
of cases where indirect methods outperformed direct approaches. Therefore, we cannot simply discard 
the indirect approaches in favor of the direct approaches. 
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APPENDIX 1 

There are numerous considerations to take into account when selecting an objective function param- 
eterization for a direct method. In this work we break down the types of parameterizations into two 
categories called the Feasible Iterate Approach, and the Infeasible Iterate Approach. In the Feasible 
Iterate Approach the cost function is parameterized in such a way that the rendezvous constraints are 
satisfied implicitly for each cost function evaluation. Using a Feasible Iterate Approach allows one to 
choose between both constrained and unconstrained numerical optimization packages. However, we must 
provide a robust way to solve the TPBVP. In the Infeasible Iterate Approach the rendezvous constraints 
are not necessarily satisfied for each cost function evaluation. The rendezvous conditions are satisfied 
upon convergence of the numerical optimization routine. The strength of the Infeasible Itera.te Approach 
is that we do not have to provide a robust TPBVP algorithm. However, it is necessary that we use only 
constrained optimization packages. Hence, in the Feasible Iterate Approach we are in a sense applying 
a change of variables to convert a constrained problem into an unconstrained problem. In this appendix 
we present one possible parameterization for each category to illustrate some differences in the methods. 
One possible parameterization of the Infeasible Iterate Approach is 

Given: (r 0 ,v 0 ,£ 0 ) and (46) 

Choose the independent variables: 

ti 


t = l,2...n (47) 

j — 1,2, 3... .n (48) 
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(49) 

(50.1 


r*i - f (r 0 . v 0 . t 0 . ti) = 
r „ - f(r/. v/.//. f„) = 

Given Eqs. (47) and (48). the entire maneuver sequence is defined. By solving Lambert's problem for 
each trajectory segment we can calculate the total AV. 

Although many parameterizations of the Feasible Iterate Approach are possible, we only present one 


here. One possible parameterization is 

Given: (r £> . v 0 . i 0 ) and (r f . v r . t j ) l 51 1 

Choose the independent variables: 

U i = 1,2. .... n (52) 

AVj j = 1,2. ..~*n - 2 (53) 


For the this parameterization the entire maneuver sequence is determined and we can solve for the total 
AV’ and satisfy the boundary conditions simultaneously. We first determine iq and vj" from Eqs. (13) 
and (14). The next step is to solve n - 2 initial value problems by iterating on the following algorithm 

for i = 1 to n - 2 


v, + 

= V-+AV) 

(54) 

r i+l 

= f(r i,v?,ti,t i+ i) 

■ (55) 

V i +3 

= g(r„v+, *i,t i+1 ) 

(56) 


end 

Finally we ensure the rendezvous conditions are satisfied by solving Lambert's problem for r n _i, r n . and 
At = t n — tn—i ■ V^ith the solution for Lambert’s problem we can solve for and Ai/„ and we then 

solve Eq. (19) for the total AV . There are many more possible choices for independent variables that 
fall under the Feasible Iterate Approach. Presenting all of the possibilities is beyond the scope of this 
work. However, it is worth mentioning that other Feasible Iterate Approaches are likely to be a hybrid 
of the method described by Eqs. (20), (21) and (22) or the method described by Eqs.(51), (52) and (53). 

Choosing a specific parameterization for the objective function is nontrivial. Each of the parame- 
terizations discussed above have some strengths and some weaknesses. Although a detailed comparison 
of the Feasible Iterate Approach and the Infeasible Iterate Approach is beyond the scope of this work, 
preliminary comparisons suggest that the Feasible Iterate Approach performs better. This is expected 
because it is often better to convert a constrained problem to an equivalent unconstrained problem if an 
appropriate change of variables is possible. Therefore we have chosen to consider only Feasible Iterate 
parameterizations. Choosing a specific parameterization of the Feasible Iterate Approach is also non- 
trivial. For this paper we choose to parameterize the objective function using the independent variables 
given in Eqs. (20), (21) and (22). The justification for choosing this parameterization is discussed in a 
previous section. 
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